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Abstract 

We present an algorithm for solving the radiative transfer problem on massively parallel computers using adaptive 
mesh refinement and domain decomposition. The solver is based on the method of characteristics which requires an 
adaptive raytracer that integrates the equation of radiative transfer. The radiation field is split into local and global 
components which are handled separately to overcome the non-locality problem. The solver is implemented in the 
framework of the magneto-hydrodynamics code FLASH and is coupled by an operator splitting step. The goal is the 
study of radiation in the context of star formation simulations with a focus on early disc formation and evolution. 
This requires a proper treatment of radiation physics that covers both the optically thin as well as the optically thick 
regimes and the transition region in particular. We successfully show the accuracy and feasibility of our method in a 
series of standard radiative transfer problems and two 3D collapse simulations resembling the early stages of protostar 
and disc formation. 
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1. Introduction 


Radiative feedback plays a crucial role in the process of star and disc formation, the evolution of circumstellar 
discs and the thermodynamics of the interstellar medium (ISM). Massive stars emit a large number of energetic UV 
photons and strongly d etermine the structure of giant molecular clouds (GMCs) by creating large bubbles of ionized 
gas (HII regions) (e.g. Peters et ^ 20ld: Walch et al. . 2012 : Dale et al. . 2013h . On smaller scales, low mass and 
intermediate mass stars also significantly influence their surroundings by radiative heating. By increasing the frag¬ 
mentation scale, radiative heating can completely inhibit fur t her fragmentation in a radius of several AU and prevent, 
e.g., the formation of a binary system ( Price and Bate . 20101) . Offner et al. (2009) investigate the initial mass function 
(IMF) and the star formation rate (SFR) by comparing 3D hydrodynamical simulations of low mass star formation 
with and without the effects of radiative transfer. They find that the thermal support of a protostar’s accretion luminos¬ 
ity suppresses further fragmentation in the cloud core as well as in the protostellar disc. The SFR in their simulations 
is about half the value of the simulations wit hout radiativ e transfer and the mass distribution of protostars of very low 
mass (M* < 0 .IMq) is significantly reduced. iBatel (120091) finds similar effects. 

Regarding the formation and evolution of circumstellar discs , radiative feedback is indispe nsable to understand their 
fragmentation behaviour, thermodynamics, and morphology ( Chiang and Goldreich . 19971) and to m odel the infrared 
excess observed in their spectral energy distributions (SEDs) (e.g. Dullemond and Monniei . 2010l) . The initial for¬ 
mation of massive discs during th e Class 0 phase has been investiga t ed using hydrodynam i cal and magnetohydrody- 
namical (MHD) sirn ulati ons (e.g. Torkeetak, 1993^ Mellon and Lil 2008 : Machida et all 2010l : Peters et ^ 2010l : 
Seifried et al. , 2011 ), and Seifried et alJ ( 2013h emphasize the importance of turbulence to explain the formation of 
Keplerian discs even if strong magnetic fields are present. 

Despite a large number of studies, the actual transition from the early self-gravitati ng protostel l ar dis c (Class 0) to 
the Keplerian protostellar disc is still poorly understood. Recent observations (e.g. iTobin et al.L 1201 2h indicate that 
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Kepler ian discs might form very early during the protostellar evolution and the analytic study by iForgan and Rice 
( 2013 emphasizes the effects of radiative processes. However, the effects o f radiative transfer have usually been ne 


glected in MHD simulations so far or wer e substantially approximated (e.g. Yorke et ah . 1993t Mellon and Lil 2008 


Machida et al.i 20ld Seifried et ah . 2011 ). The self-consistent modelling of the formation and early evolution of pro¬ 


tostars and protostellar discs therefore creates the need for numerical methods to make 3D radiation MHD simulations 
feasible. 

In this context, radiative transfer is a rather costly computation compared to solving Euler’s equations. The reason for 
this is that the timescale of radiative transfer is usually much shorter than those of hydrodynamics and MHD because 
of the large speed of light compared to the sound speed of the gas in, e.g., a molecular cloud or the characteristic 
Alfven wave speeds of the magnetic field. The short timescale on which radiation emerges throughout the complete 
computational domain makes radiative transfer a highly non-local problem compared to MHD which is determined 
completely by local thermodynamic properties of the gas. In this sense, hydrodynamics and radiative transfer are 
two very different numerical tasks and very challenging to solve consistently. Modern Eulerian MHD codes like 
EL ASH ( Ervxell et al.L 2000l) mostly solve the Euler equations on a grid with adaptive mesh refinement (AMR) to 
resolve fluid features on a wide range of length scales. These codes are parallelized by subdividing the computa¬ 
tional domain into several subdomains each containing a fixed number of cells. Since the Euler equations describe 
local fluxes of mass, momentum and energy, all subdomains can be handled in parallel during a hydrodynamical time 
step. Between the time steps, boundary values of the subdomains are exchanged using the Message Passing Inter¬ 
face (MPI) for communication. In contrast, characteristics based radiative transfer codes are usually designed very 
differently. Instead of domain decomposition, these codes are parallelized exploiting the formal independence of the 
radiative transfer equation (RTE) on the solid angle. Resolving the anisotropy of the radiation field accurately requires 
a large set of characteristics each covering a discrete opening angle of the An unit sphere. If all radiative quantities 
are assumed fixed during one solution step, characteristics of different directions can be computed independently of 
each other which makes it ideal for parallelization. However, the spatial information of the computational domain 
with all radiative quantities has to be available to the each processor computing a certain number of characteristics 
on the solid angle grid. This can be a severe drawback in terms of memory requirement if high spatial resolution is 
required or a large number of frequencies or both (e.g, synthetic stellar spectra). Solving both Euler’s equations and 
the RTE consistently requires careful approximati ons to the radiat i ve tran sfer problem to make the coupling of an 
MHD code with a radiative transfer code feasible, van Noort et al. ( 2002h present a radiation solver that is coupled 
to a hydrodynamical code using AMR and domain decomposition in 2D. The radiation solver uses short character¬ 
istics (SC) for integrating the RTE while boundary values are communicated between Lambda iteration steps. The 
focus of this approach lies on modelling the dynamics of scatteri ng dominated stellar atmospheres. The SC ap proach 
allows for a fast converging Gauss-Seidel iteration scheme (e.g.. iTruiillo Bueno and Eabiani BendichcJ,ll995b . while 
non-local contributions have to be communicated by a suc cessive exchange of boundary values b etween subdomains. 
This approach was also extended for 3D simulations (e.g., Havek et^D. 201(i Davis et al. . 2012 ). However, while the 
Gauss-Seidel short characteristics approach is well suited for highly scattering domi nated regimes, it introduces a lot 
of numerical diffusion because a large number of upwind interpolations is necessary. iRazoumov and CardaUl (120051) 
implement a method that is as computationally cheap as the SC method hut less dilTusive. They create rays 
on each refin e ment level separately while their approach is fully threaded hut not MP I parallelized. Recently, 


Tanaka et alJ (120141) parallelized this approach using the multiple wavefront method hy iNakamoto et alJ (120011) 


based on a carefully chosen calculation sequence on a spatially decomposed domain. This method requires 
successive communica tion of boundary values. A similar approach is used with long characteristics (EC) in 3D by 
Heinemann et al.i ( 20061) without AMR. 


Another approach for including radiative transfer in hydrodynamical simulations is based on the moment 
equations (the angular integrated RTE) of the zeroth, first and second moment of the specific intensity. A 
moment-based scheme does not necessarily require to integrate along large sets of characteristics, however, the 
anisotropy of the radiation field has to approximated reasonably in order to close the set of moment equations 
for the mean intensity, r adiative fiux and pres sure. A possible closure relation is the Ml-closure used, e.g., 
in the HERACLES code dConzalez et all l2007h . The closure relation can also be explicitly calculated using, 
e.g., a characteris tics based approach which is known as the Variable Eddington Tensor (VET) method (e.g. 


.liang et alJ. 120121) . A common approach for star formation simulations is the diffusion approximation of the an¬ 


gular moment equations which assumes the radiation field to be completely isotropic. In regions of high opacities 
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X, the diffusion approximation is an expansion of the specific intensity in which all terms cx 1 are neglected in 
the RTE. This leads to Eddington’s approximation in which the isotropic radiation pressure is proportional to the 
radiation energy density. The moment equations of the radiative intensity themselves then form a set of hyperbolic 
equations, like Euler’s equations. However, since those two hyperbolic systems would still have to be handled on their 
individual timescales, one can even make one further step and neglect the time dependence of the radiation flux by 
assuming it to be proportional to the gradient of the radiation energy (Eick’s law). The moment equations can then 
be combined into a single diffusion equation for the energy of the radiation field. Because the flux in the diffusion 
approximation lost its finite propagation speed, one has to introduce a flux-limiter to a void unphysical propagation 
speed s depending on the actual opacity. This, flux-limited dijfusion approximation (ELD) ( Levermore and Pomranina 


19811) has been successfully used in radiation hydrodynamical star formation simulations coupled within Eulerian grid 


codes (e.s.lStone et al.. 19921 Krumholz et al.. 20071: Commercon et al.. 2011: Flock et al.l 

20131 Zhang et al. , 2013; 

Brvan et al.l 20141) as well as smoothed particle hvdrodvnamics (SPH) codes fe.g. Bate et al. 

2013). However, the dif- 


fusion approximation is only valid in optically thick regions where the radiation field becomes isotropic. iKuiner et al 


( 2010h have shown the significant drawbacks of exclusively using ELD in the transition regions from optically thick 
to optically thin regimes where the radiation field becomes highly anisotropic. Recent efforts have bee n made to 
combine raytracing methods with ELD solvers ( Kuiper et al.l 2010t Elock et al. , 2013 ; Klassen et al.l 20141) to handle, 
at least, primary stellar or protostellar radiation separately from the ELD approximation and to avoid the stellar flux 
from diffusing into shadow regions. 

Einally, Monte Carlo (MC) methods have become increasingly popular during the last decade, especially in post¬ 
processing MHD simulations. The MC method is a statistical approach and treats individual photons or photon pack¬ 
ages by following its propagation path and comput ing absorptio n, emission and scatteri ng probabilities. Several ad- 
vances have been introduced, e.g., photon peel- off (Lucy, 1999 ). immediate reemission ( Biorkman and Woodl 200 ih 


and diffusion approximations (IMin et al.L l2009h which make the MC method a powerful tool to calculate synthetic 


spectra, SEDs or polarization maps from the outcome of MHD simulations. The angular and frequency resolution are, 
in principle, unlimited since the direction of propagation of a photon package and its frequency are chosen randomly 
from a continuous probability function. In that sense, the MC method always gives a quite reasonable result even in 
the limit of a small number of photon packages while a low resolution shows mainly up as statistical noise in the solu¬ 
tion. But the statistical approach also has a severe drawback since we do not know in advance the exact path a photon 
package will travel, and how and when it is emitted or absorbed. Therefore, it is extremely difficult to implement on a 
decomposed domain. MC methods are extremely successful in post-processing the outcome of MHD simulations but 
are r ^ely used in combina tion with hydrodynamical simulations. Those approaches which does include MC methods 
(e.g. Acreman et al.L 2010h are fairly restricted in their spatial resolution of the AMR grid, since each processor has to 
get a copy of the complete computational domain to be able to follow the path of an arbitrary photon package. Eor our 
approach, we therefore choose a discrete ordinate method using characteristics to integrate the RTE which requires a 
raytracer that works on an AMR grid with domain decomposition. 

The paper is organized as follows: In Section |2] we give a brief introduction into the theory of radiative transfer as 
far as it concerns our method and we describe in detail the method of hybrid characteristics. We also describe the 
coupling between our radiation solver and the ELASH code. In Section [3] we show results from test calculations we 
performed to investigate the accuracy of our radiation code as well as the coupling to the ELASH code. In Section|4l 
we present results from 3D radiation hydrodynamical collapse simulations and the parallel scaling performance of our 
code is described in Section |5] In Section |6] we discuss our results and put it into context with other state-of-the-art 
radiation transfer methods. 


2. Theory and Numerics 

In this Section, we describe the theory of radiative transfer (RT) that forms the basis of our solution method as 
well as the numerics. We describe the hydrodynamics only as it becomes important in the coupling with the ra diation 
solver. Eor a more detailed description of the ELASH code and its capabilities we refer to Ervxell et al. ( 2000l) . 


2.7. The Equation of Radiative Transfer 

The theory of radiative transfer in this section is based on the work by Mihalas and Weibel Mihala^ ( 1984 ) in 
the limit of geometrical optics. The energy of the radiation field is described by a scalar field of specific intensities 
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I(x, n(0, (p), v), where x is the position in space, 6 and cp define the direction of propagation n, and v is the frequency. 
The radiative transfer equation (RTE) describes the change of the specific intensity during its propagation in a medium 
which is determined by an energy balance between emission and absorption processes. It reads 


+ n • V/y — T]y ^y/y, 

c ot 


( 1 ) 


where rjy denotes the emissivity (energy volume density per unit time and solig angle), is the extinction coefficient 
(1/length) and c is the speed of light. The specific intensity denotes the radiative energy flux per solid angle dO. - 
sin6d6d(p, thus in vacuum, it is constant along a line of sight. Interaction processes between radiation and matter 
determine the extinction coefficient ;^fv. However, for this work, we solve the time independent RTE since we assume 
the radiation field to emerge instantaneously throughout the entire computational domain during a hydrodynamical 
time step. Eurthermore, we use the definition of the source function Sy - rfyjxy and rewrite the RTE in terms of the 
optical depth without explicit frequency dependence: 


dl{n) 

dT{n) 


^S- /(n). 


( 2 ) 


This form of the RTE describes the propgation of the specific intensity along a specific line element ds in the direction 
n and the optical depth element dr - xds respectively. This requires a proper paramerization depending on the 
coordinate system and, hence, the definition of n ■ V. Note that the RTE in the form of Equation (|2ll becomes a ID 
ordinary differential equation. However, for the numerical solution in 3D, the optical depth element dr is discretised 
and parameterised in Cartesian coordinates and the solution is obtained by integrating the RTE on a solid angle grid. 
The source function 5 is a more general form of Kirchoff’s law. It describes the ratio of emission and extinction 
of radiative energy and allows arbitrary contributions from thermal emission as well as scattering contributions. In 
fact, the complexity of the model and the solution of the RTE depend strongly on the source function we choose to 
accurately describe the current radiation transfer problem (e.g. local thermodynamic equilibrium (LTE), non-LTE 
(NLTE), grey or non-grey, anisotropy, dust continuum radiation, line transfer, etc.). Describing the complete radiation 
field would require 6 dimensions, which makes it an extremely challenging task to compute and store a 3D solution. 
In order to handle the radiation field numerically, the intensity is only computed on the fly and accumulated in form 
of the solid angle averaged mean intensity 

IdQ. (3) 

47r J4;r 

The mean intensity is the zeroth moment of the specific intensity and closely related to the radiative energy density 
Ey - Depending on the model setup, the source function usually depends on the radiation field itself and the 
mean intensity becomes a part of the source function, which makes the RTE an integro-differential equation. In order 
to find a self-consistent solution, one has to invoke an iteration scheme of some form. Eormally however, the RTE 
from Equation (|2]i can be solved by the formal solution: 


/(T2) = /(Tl) 


H- 



S (t) ^^dr. 


(4) 


The formal solution describes the intensity propagation along a line element with the optical depth At = T 2 - ti . It 
contains the incoming intensity /(ti), which is partially extinct and additional energy from emission processes. The 
RTE and the formal solution are linear in the intensity and allow us to split the radiation field in as many components as 
the solution method requires. This is a crucial part in our approach of solving the RTE on a decomposed computational 
domain such as the adaptive grid embedded in the ELASH code. 


2.2. Numerical Radiative Transfer 

Integrating the RTE along a set of rays of different directions n using Equati on dUl is ba s ed on the method of 


characteristics. It was first introduced into the radiation transfer community by Olson et al. ( 19861) . The RTE is 


integrated for each cell in the computational domain and each direction by computing a stepwise formal solution 
along a ray, or long characteristic (LC), according to the discretised formal solution 


h = li-x exp(-AT,_i) H- Mi, 


(5) 
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Figure 1: The staggered RT-grid (solid lines) defined by the 
cell centers of the underlying finite-volume hydro-grid (dashed 
lines); a long characteristic at an arbitrary direction is shown, 
which integrates the RTF for the hydro cell-center at the upper 
right corner of the domain. 


where At,- is the hnite optical depth element given by a piecewise linear interpolation 

= ^0,-1 

Xi is the opacity at the discretization point i, on the characteristic. A/, is the discretised counterpart to the integral in 
the formal solution (Equation (|4]i) and is solved by either a linear or parabolic interpolation according to 


+ PiS i -t y/iSz+i. 


(6) 


The co efficients a,-, /?, and y, depend on the optical depths between s,_i, s, and s,+i. They are given in lOlson et al 
( 19861) . Figure [T] shows the geometrical situation of a characteristic passing through a homogenous grid at an arbi¬ 
trary direction ny(0,0). Since the opacity and source function are stored in the cell centers of the finite volume 
FLASH grid Mashed lines), they are assumed to he constant inside the cell. However, since we use a parabolic 
interpolation 0 of the source function integral, we introduce a point-hased RT grid which is based on the cell 
centers of the FLASH grid. These cell centers define the vertices from which we interpolate the values at the 
intersection points of the ray bilinearly. Consequently, the point-based RT grid is staggered by half a grid cell 
since the ray does not intersect with the grid faces of the finite volume grid but with the faces of the point based 
grid defined by the FL ASH cell centers. The chara cteristic is traced on the RT-grid using the fast voxel traversal 
algorithm introduced by Amanatides and \\^ (1987). The opacity and the source function 5, at the intersection 
points of the characteristic with the RT-grid are interpolated bilinearly from the adjacent vertices. 

While the RTE is integrated along each direction n, the mean intensity is computed by accumulating all intensities: 


•' = 3; 


(7) 


simple cell based linear interpolation leads to a loss of radiative energy in optically thick regimes, especially if iteration is required (e.g. 
Dullemond 2013, http://www. ita.uni-heidelberg.de/~dullemond/lectures/radtraiis_2013/index. shtml Chapter 3.8.4 and 4.4.8 
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which requires a discretization of the solid angle Q. If no information about the anisotropy of the radiation field is 
available, one should choose a homogeneous discretization which is a non-trivial problem if one considers spherical 
coordinates on the unit sph ere. For this purpose , we use the HEALPix (Hierarchical Equal Area isoLatitude Pixeliza- 
tion) scheme introduced by Gorski et al.l ( 20051 ) . HEALPix ensures an optimal discretization of the unit sphere into a 
number of finite solid angles AQ. The discretization is based on 12 base pixels which are subdivided depending on 
the required resolution level. Consequently, typical numbers of directions for the integration of the specific intensity 
are Apix = 12,48,192 or 768 ( AppendixB| l 

Characteristics based radiative transfer is the attempt to approximate the radiative interaction of each cell with each 
other cell in the computational domain. Although the method of long characteristics is very accurate in doing this, it 
is rather inefficient as it requires to shoot a large number of rays for each cell to sample the radiation field accurately 
in 3D. An alternative is to use a short characteristics (SC) approach, in which only neighbouring cells are used to 
interpolate incoming intensities from different directions. This requires to sweep the cells in an ordered fashion to 
make sure that all intensities which are required for interpolation are available. The SC approach introduces numerical 
diffusion because of the large number of interpolations involved but reduces the cost of the RT calculations by a factor 
of ric (the total number of cells involved). Either way, the method of characteristics must invoke a raytracer, which 
samples radiative interactions between arbitrary regions in the computational domain. 


2.3. Ray tracing on the decomposed AMR Grid 

The parallel design of the ELASH framework, in which our solver is currently implemented, forbids to trace rays 


over th e entire domain as it is necessary for the method of characteristics. ELASH invokes PARAMESH (lOlson et al 


digQgt) ). and lately also the CHOMBO libr arjH, for implementing an adaptive mesh refinement (AMR) grid. Paramesh 
uses a block structured AMR mesh, in which the fundamental data structure is a block containing cells which are logi¬ 
cally indexed by a coordinate triple (i,j,k). The entire computational domain consists of a number of blocks of different 
physical sizes ordered hierarchically in an octree data structure. Blocks at the bottom of the tree structure, called leaf 
blocks, contain valid data and they cover the entire physical size of the computational domain. PLASH allows for 
massively parallel computation by invoking the Message Passing Interface (MPI) for the communication of ghost cell 
information between the blocks. Optimal load balancing is guaranteed by splitting the AMR tree equally between all 
available MPI tasks to ensure that each task receives more or less the same number of leaf blocks. E.g., the AMR 
tree of a star formation simulation typically requires more than 10 levels of spatial resolution with up to several 10^ 
blocks each containing 8^ cells, which is only made possible (in terms of cpu-time and memory requirements) by the 
parallelization described above. 

The method of characteristics stays in direct contrast to the spatial parallelization of the AMR grid (Pigure|3. How¬ 
ever, in order to account for non-lo cal coupling of th e radiation field, w e adapt a raytracing technique originally 
developed by Riikhorst et al.l ( 2006h and improved by Peters et d1 ( 2010l) . which uses a combination of local long 
characteristics and a global ”short-characteristics-like” interpolation of outgoing intensities from the decomposed 
domains of the AMR grid. The basic idea is to split the radiation field in two components: 

• LA local component which uses long characteristics to compute only radiative contributions to both the cells 
inside a block as well as contributions that leave the block {face values). The computation is done in parallel 
and in accordance with the design of the block structured AMR grid. 

• 2. A global component which is computed by communicating and accumulating face values (see EigureO. 
This step invokes raytracing over the block structure of the AMR tree and a linear interpolation of face values 
very similar to the SC method (but on the level of subdomains). After the communication of face values and 
the tree hierarchy, this step is also done in parallel. 

This approach, called hybrid characteristics, only needs to communicate the face values of the blocks and information 
about the AMR tree hierarchy but no 3D data. By this, the amo unt of communicated data is reduced significantly. 
Originally, this method was developed by Riikhorst et al. ( 2006h to compute column densities only with respect to 
point sources for UV ionization. The original method requires to communicate the whole AMR tree structure at the 


■ https://seesar.Ibl.gov/anag/chombo/ 
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Figure 2: Example for a 2D AMR grid distributed over two processors without shared memory. The bold rectangle shows 
the boundary of the whole computational domain. Thin lines show the leaf blocks at different refinment levels that make up 
the whole subdomain a processor is assigned to. Raytracing through the domain is obviously restricted to the subdomain 
each processor is assigned to. 


highest level of spatial resolution during the raytracing step on the AMR block structure. This stands in contrast to 
the parallel design of the FLASH code and re stricts the available range of refinement levels of the AMR tree substan¬ 
tially because of the large memory overhead. Peters et H] ( 2010ll add some major improvements to the algorithm by 
introducing a walk through the AMR tree, which only requires the communication of basic AMR information and 
conserves the idea of shared memory parallelization. 

However, the method was, originally, restricted to compute column densities along rays which originate at a 
certain point in the grid and used to represent, e.g., a stellar source. For this work, we removed this restriction 
and implemented a radiative transfer framework which is able to compute the radiation field independently 
from any point source hy solving the RTF for large sets of characteristics along parallel rays. By comhining 
our improve ments with the or iginal meth od, the solver can n ot only account for the primary emission hy point 
sources (as in iRiikhorst et alJ (200^ and IPeters et ^ dlOlOb ) hut also for the reemitted, diffuse component of 
the radiation held. Figure |4] shows a 2D example of a simple test setup with an irradiated central density clump 
using AMR. From the figure we can see the ability of the method to create sharp shadows and to transport incoming 
radiation over the entire domain. 


2.4. Coupling to the FLASH Code 

Since our method is implemented in the FLASH framework, it is straightforward to couple the radiative transfer 
module to the hydrodynamical and MHD modules of the FLASH code. The coupling is done by accounting for 
radiative emission and absorption processes, which are determined by the thermal emission opacity = Kep and the 
thermal absorption opacity = Kap. The opacities are calculated from mass specific cross sections Ke and Kg. Note 
that the total extinction coefficient;^, which is used for the solution of the RTF, may include an additional scattering 
opacity and therefore = = Yn+Ys:. The coupling of both the ra diation and the MHD solver is achieved by computing 
a source term according to iMihalas and Weibel MihalasI (119841) which describes the total net gain or loss of energy 
due to radiative heating and cooling. It reads 


Qrad = -S) ^ AnXaiJ - B). (8) 

This source term is computed from the time-independent solution of the radiation field as described in the previous 
sections and it is coupled to the MHD integrator in an operator splitting step. Hence, the set of compressible MHD 
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(c) interpolation 




Figure 3: The basic steps of the hybrid characteristics method for parallel rays (compare to lRiikhorst et alJ j200d) ) in a 2D 
AMR domain that is refined from left to right. Bold lines show the boundaries of the patches at dilferent AMR levels (in 
FLASH, these patches are called blocks and are distrlbted equally over the available number of MPI tasks). In this example, 
each block contains 4x4 cells (indicated by thin lines), a) Local contributions as calculated with long characteristics, b) The 
outgoing face values which are communicated. In fact, we communicate all face values even though they might be part 
of the same subdomain of a certain MPI process since we need them in the following interpolation step, c) Example for 
the linear interpolation of face values for a particular target cell after the communication step. The linear interpolation 
requires to weight the face values from two rays. The weights depend on where a certain ray segment starts at the boundary 
of a block (which is a 4x4 cell patch at a certain refinement level) or subdomain respectively. 
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Figure 4: The specific intensity and the optical depth computed diagonally in the xy-plane with a central density clump. The source 
function is set to unity at the left and bottom outermost boundaries and zero everywhere else. The opacity of the central clump is 
one order of magnitude higher than the ambient opacity. The grid indicates the block structure of the AMR grid, units are arbitrary. 


equations in dimensionless form including gravitation and radiative energy exchange are those of 


continuity 
dp 


dt 


+ V ■ (pv) = 0, 


momentum conservation 
d(p\) 


dt 


+ V ■ (p\ ® V + p,1 - B ® B) = pg 


energy conservation 
dE 


dt 


+ V ■ (\{E + p.) - B(V ■ B)) = pv ■ g + 2rad, 


and the induction equation 

5 - V X (V X B) = 0, 
dt 


(9) 


( 10 ) 


( 11 ) 


( 12 ) 


with the gas velocity field v, the magnetic field vector B and the gravitational acceleration g. p» is the total pressure 
and E the total energy density of a fluid element containing magnetic contributions according to 


P*^P+—^ 
1 




E = -pu +emp+ 


(13) 


(14) 


with the gas density p, the thermal pressure p and the internal specific energy ei„f 1 denotes the unity matrix. The 
MHD equations are closed by an ideal gas equation of state (EOS) which relates the internal energy to the thermal gas 
pressure according to 

p^ir-l)pemf (15) 

We assume y = 5/3 which corresponds to a mono atomic (hydrogen) gas. The temperature is also related to the 
internal energy by: 

(16) 

pmp 
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where kt is the Boltzmann constant, irip the proton mass and yu the mean molecular weight of the gas. 

Note that we solve the equations of MHD and RT successively by an operator splitting step and not simultaneously. 
Furthermore, for the following test cases, the thermal pressure dominates the hydrodynamics and it is several orders 
of magnitude larger than the radiation pressure, which we therefore neglect in the momentum equation (fTOl) . However, 
since our method explicitly computes the angular dependency of the radiation field, it is straightforward to couple it 
into the MHD equations. 

2.4.1. Choosing the Time Step 

The current coupling is done by an update of the internal gas energy eint and temperature T respectively. Since we 
solve the time independent RTF, there is no update of the radiative energy or the source function during the solution 
of Euler’s equations. Instead this is done in the following time step when the gas quantities have been updated. The 
update of the internal energy is done explicitly by 


Aeint = Af 2rad- (17) 

Due to the explicit update, we have to make some restrictions on the time step. The radiation field does not have 
an explicit influence on the CEL time step since the energy update is done after the solution of the MHD equations. 
Instead, we compute a cooling time step which is chosen if it is shorter than any other time step from a FLASH 
module. The cooling time step is chosen so that the energy contribution Aeint does not exceed a fixed percentage of 
the internal energy. Otherwise, if the time step Af is chosen too large, the total radiative energy could become negative 
(e.g., Ae^ > e). This leads to the following time step restriction: 


Afc = min(--^ll^^] kc AfcFL, ( 18 ) 

where kc determines how much change in the internal energy is allowed, Afcpu is the CFL time step, and i 
denotes the indices of all grid cells in the computational domain. Because of the explicit energy update, the 
cooling time step is usually shorter than the CFL time step. So far, there is no suhcycling involved and the 
FLASH code chooses a global minimum time step from all physics modules involved (including, e.g., self¬ 
gravity). The cooling time step highly depends on the absorption coefficient ;t' since it determines the optical depth 
of the medium and how much radiation is absorbed and emitted during a single time step. Typically the choice of 
0.2 > kc> 0.01 is convenient as it produces accurate results (Section [331 l and time steps about one oder of magnitude 
lower than the CFL time step. 


2.5. The Lambda Formalism 

Computing the radiation field in the form of the mean intensity in Equation O requires a formal solution of 
the RTF in the way described above. Usually, this task is described in a rather compact form by using the Lambda 
Operator. 

7 = A[5]. (19) 


Formally, the Lambda operator for one cell in the computational domain contains the radiative contributions from each 
other cell. The construction of the operator would require to explicitly calculate the radiative coupling between a cell 
and each other cell. But this is far too costly concerning computation time and memory requirements. Instead, we do 
not construct the operator but we approximate the Lambda step from Equation ( fT9l l by using the formal solution from 
Equation (|4]i to compute the radiation field J from the source function S in the way described above. The accuracy of 
this approximation in a 2D or 3D computation depends crucially on the angular resolution, since it determines whether 
we actually ’’hit” each other cell during the angular integration of the mean intensity or not. We avoid this prohlem 
partly hy calculating the radiation from point sources (e.g. a stel lar source) explicitly fo r ea ch cell hy com hining 
our method with the original hybrid characteristics method bv iRiikhorst et alJ ( jlOOtil) and fPeters et alJ dlOlol) . 
However, the Lambda step from Equation [19] requires that we know the source function in advance. If we take the 
temperature from FLASH’S hydro solver, we can compute the source function simply as being S - B{T) then solve for 
the radiation field, couple it back to the hydro solver and we are done. This approach assumes the gas to be in a state 
of thermodynamical equilibrium but this is, of course, not always the case. If the radiation field is decoupled from the 
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gas temperature, we do not know the source function in advance. The solution then requires some kind of iterative 
procedure to account for the non-local coupling of the radiation field with the gas. In the theory of radiative transfer, 
this iterative method is called Lambda iteration, which requires iterating over Equation (fT9]) until a self-consistent 
solution for J(S) is found. Strictly speaking, even in the LTE case with S - B(T), we have to iterate to find a 
temperature that is consistent with the internal energy of the gas since this determines the thermal emission. However, 
the Lambda iteration may need several hundreds of iteration steps, which is too costly and ineffective to be employed 
in a hydrodynamical simulation. One way of resolving this problem, is to partly solve Equation ( fT9b analytically by 
splitting the Lambda operator. These approaches, called Accel erated Lambda iteration (ALI), have been in vestigated 
and used extensively in the stellar atmosphere community (e.g. Trujillo Bueno and Eabiani Bendicho[ll995h . We have 
implemented the most simple form of ALI, the local lambda operator, to solve radiative transfer problems even in 
regions of high optical depths and strong decoupling where the classical Lambda iteration usually fails (| Appendix A|i. 


3. Tests 

In this section, we show test results from the implementation of our radiation solver. The tests include time 
independent (Sections l3.1l and l3T2l i as well as dynamical tests (Section UAT i in ID and 3D. We also show results from 
the combined LLASH/RT code in a series of ID radiative shock calculations in SectionjTT] 


3.1. ID Atmosphere 

In the first test, we compute the radiation field in a grey, isothermal, scattering dominated ID atmosphere. This 
test is typically used to verify a radiation solver’s iterative performance and accuracy in a non local thermodynamic 
equilibrium (NLTE) situation on a wide range of optical depths. It is also particularly useful to ensure that the solver 
accurately reproduces the diffusion limit in an optically thick regime, e.g., in the lower parts of the atmosphere. This 
test also requires the ALI method since the classical Lambda iteration fails to reproduce the solution in the case of 
strong scattering contributions. The amount of scattered radiation is quantified by the ratio of the thermal absorption 
coefficient to the total extinction coefficient according to 


e - 


Xa 

Xa 


( 20 ) 


where we neglect the frequency dependence and e is the photon destruction probability. The grey source function in 
the atmosphere contains a thermal part and a scattering contribution, and it reads 


1 ^ 

X Xa+Xs Xa+Xs' 


( 21 ) 


= (1 - e)J + eB, 


( 22 ) 


where we defined J - rjslxs, the thermal emission is B = rjelXa, and ps and rje denote the scattering and the 
thermal emissivity respectively. Since the atmosphere is isothermal, we assume that we know the temperature and 
normalize it so that B - 1. The crucial part in this test is to find the source function which has to be consistent with 
the mean intensity J which is 

J = i(/. + h), (23) 


where I and /+ are the down and upward (2 stream solution) integrated specific intensity respectively. Since we 
assume a uniform mass specific opacity k and constant temperature T, the intensity is only a function of optical depth 
dr - xdz, thermal emission B and the photon destruction probability e. The mean intensity is then given by the 
analytic solution 


7 = B 



exp(-Vir) ) 

1+Vi /■ 


(24) 


The density p of the model increases exponentially with distance from the upper boundary and we assume that x P 
but with e being constant. There is no incoming radiation at the upper boundary of the atmosphere at t = 0 while 
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non-LTE plane parallel atmosphere: two stream solutions 



optical depth x 

Figure 5: Scattering dominated ID atmosphere problem. The solutions from the radiation solver (symbols) are compared to the 
analytic solutions (lines) for five different photon destruction probabilities. 


at the lower boundary the incoming radiation is I - B. The resulting model atmosphere provides an exponentially 
varying optical depth t which resolves the transition region from the optically thick inner LTE-regions to the optically 
thin NLTE-regions at the outer boundary. We test the solver for a wide range of photon destruction probabilities 
from e - 10 ' to 10“^. The domain consists of 8 subdomains each containing 8 cells which results in a total spatial 
resolution of 64 cells. Figure |5] shows the results. In the outer optically thin parts of the atmosphere, the scattering 
contribution in the source function becomes dominant since radiation leaves the atmosphere. The numerical solution 
is in excellent agreement with the analytic solution. 


3.2. Hydrostatic Protostellar Disc 

Cosmic dust is one of the most important constituents of the ISM. By mass, it makes up only a small fraction 
of typically about 1%, but dust has important radiative and chemical properties. Dust particles have strong contin¬ 
uum opacities which are highly frequency-dependent. Especially in the optical regime, dust absorbs light much more 
efficiently than in the infrared regime. That is why young protostars, which are surrounded by gaseous and dusty en¬ 
velopes, are difficult to observe in the visible wavelengths but require infrared observations. Thermal absorption and 
reemission of radiation by dust (a process called reprocession) strongly determines the thermodynamical properties 
of a protostellar disc, especially in those regions where the disc is opaque to direct stellar radiation and dominated by 
thermal reemission of dust molecules. This is mainly the case near the equator of the disc because radial optical depths 
with respect to the central star are typically much larger than unity (t» » 1). Therefore, modelling the temperature 
structure requires diffuse radiative transfer to be taken into account. 
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In this test setup, we combine emission from a point source with the solution of the RTE. The goal is to determine 
the self-consis t ent te mperature structure of the gas in a protoste llar disc. The setup is based on the benchmark by 
Pascucci et al. (2004), which is based on the theoretical work by Chiang and Goldreich ( 1997h . We co mpare our so¬ 
lution s from a 3D calculation with the results from the Monte Carlo radiative transfer code RADMC-3D ( Dullemo^i^ 

2OI2I) . 


3.2.1. Thermal Radiative Transfer 

A protostellar disc combines optically thick and thin regimes, which requires the computation of primary stellar 
radiation and the thermal reemission from dust molecules in the disc. Our approach follows the i dea of splitting 
the radiation field in two components handling each separately. Following the work of Dullemon3 ( 2002h . the first 
component we compute is the extinct stellar flux. This can be handled by using the original hybrid characteristics 
method, which computes the optical depth with respect to a central stellar source (t*). The extinct stellar flux F, at a 
distance r from a star of luminosity L, is given by 


T’*(x) = 


Anr^ 


exp(-T.(x)), 


(25) 


assuming that the star can be approximated as a point source. The amount of energy per unit time that is absorbed this 
way is determined by the absorption coefficient and given by 


2(x) = T^.(x). 


(26) 


The reemitted radiation of the dust grains in the disc is treated as a secondary component of the radiation field. This 
component is computed with the general transfer algorithm using parallel rays. Assuming LTE, the dust grains will 
acquire an equilibrium temperature such that they emit exactly the same amount of energy which they absorb 




(27) 


where I is the specific intensity of the reprocessed radiation field. The first term in Equation (IZTT i accounts for the 
direct stellar radiation while the second term describes the energy of the reprocessed radiation field. The transfer 
equation for reemitted radiation by dust grains is 


dT 



n 


(28) 


Hence, the source function in this setup is the frequency integrated thermal emission from dust grains S - 242 .7^4 "pjjg 
task at hand is to find a temperature that is consistent with the coupled set of Equations (IZTT i and (1281) . This is done by 
iterating the equations until convergence is reached (Lambda-iteration). 


3.2.2. The Disc Model 
For the simulation setui 


we ar e following the benchmark test of iPascucci et al.l (120041) which resembles flared 


;tup we an 
i^ ll997h . 


disc dChiang and GoldreichLll997h . The idea is to define a radial gas surface density distribution and to assume that the 
vertical density structure is only determined by the hydrostatic equilibrium in the vertical direction. The gas density 
distribution is given by 

p(r,z) ^ PQ Mr) fiir). 


fiir) 

Mr) 

h{r) 





(29) 
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Figure 6: The dust density in the xz-midplane for the Pascucci benchmark for a total optical depth of Tjisc = 1. 


where r is the radial distance in the disc midplane, z is the height above the disc, and po is the gas density in the 
midplane at r = = 500 AU and z = 0. The outer disc radius is defined by rout = 1000 AU = 2 id and we crop 

the disc at an inner radius r = rin. Zd determines the height of the disc which we choose to be 0.25 r^ consistent 
with IPascucci et alJ (l2004li . We choose the central source to have solar properties with M, = 1 Mq, R* = 1 Rq and 
Tt - 5800 K. We u se a grey opacity at the visible wavelength of T = 550 nm from the opacity tables used in 
Pascucci et al. (2004) {k - 8736 cm^g"'). 


In contrast to the Pascucci benchmark, we perform our calculations in 3D instead of 2D. Therefore, we can not directly 
compare our results to the Pascucci results but instead use the results from RADMC-3D as a reference. We perform 
calculations for three cases of po so that the total radial optical depth of the disc in the midplane varies from Tdisc = 1, 
T'disc = 10 and Tdisc = 100. We do not explicitly distinguish between a gas and a dust temperature and assume both 
to be tightly coupled and the dust density is defined as a fixed fraction of the gas density (1%). The dust density 
distribution through the xz-midplane of the disc setup for the optically thin case (Tdisc = 1) is shown in Figure |6l 
The linear spatial resolution varies over 4 refinement levels from Aa: = 31.25 AU in the outer regions to Ajc = 1.953 AU 
in the center of the disc. The solid angle integration is performed using 768 directions (nSide=8). 
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3.2.3. Results 


The resulting temperature structures and averaged midplane profiles are shown in Figure |7] As it turns out, the 
accuracy of the solution is very sensible to the spatial resolution of the inner edge of the disc at r = rin which 
is a result of discretizing the inner circular rim on a Cartesian grid. Therefore, we increase the inner radius from 
rin = 10 AU, 20 AU to 40 AU for the three different setups to guarantee sufficient resolution at the point where the disc 
becomes optically thick. 

In the optically thin case (Tdisc = 1), the midplane temperature is almost entirely dominated by the direct illumination 
of the central source. In the optically thick cases, the midplane temperature is dominated by the reprocessed radiation 
from dust in the photosphere of the disc, which is directly illuminated by the central source. At the point where the 
disc becomes optically thick, a bump emerges in the temperature profile since the dust distribution becomes dense 
enough to absorb a considerable amount of radiation from the central source. Our results are in excellent agreement 
with t he reference computed by RADMC-3D and within the 10% range of the results from the different codes used 
for the lPascucci et akl (120041) benchmark. 

However, the temperature structure iu the left pauel of Figure|7]is seusitive to the augular resolutiou. Although 
the raytracer takes care of the kuowu primary stellar radiatiou, the solutiou iu the outer regious depeud ou 
whether the reprocessed radiatiou from the hot iuuer rim is accouuted for correctly. Especially iu the optically 
thick case (Tdisc = 100), siugle rays become visible iu the temperature strucutre eveu for a large uumber of 
directious (Apix = 768) siuce uot each cell is correctly couuected to the hot iuuer rim iu terms of radiative 
exchauge. A larger uumber of directious is very costly, but au alteruative is to also model the emissiou from 
such ’’hot spots” as a part of the primary emissiou. The problem is to ideutify these hot spots iu the domaiu 
siuce they cau uot be represeuted by poiut sources or siuk particles. But ouce these regious a re ideutifled, 
their emissiou cau be haudled by au iuverse raytracer similar to the approach for poiut sources dPeters et all 


2010l) . However, this approach requires au adaptive augular grid while our approach is ouly capable of usiug a 


homogeuous augular resolutiou at the momeut (AppeudixB). 


3.3. Diffusion Test 

In this section, we show results from a time dependent radiative transfer calculation. Solving the time dependent 
RTF on the timescale of the speed of light would lead to time steps far too small for the use in a hydrodynamical 
simulation on astrophysical scales. However, since we are not interested in the dynamics of the propagation of the 
radiation itself but in its contribution to the energy budget of the gas, we assume the hydrodynamical timescale to be 
much larger than the timescale on which radiation is transported. This means that the radiation field emerges instanta¬ 
neously everywhere, and we assume the solution of the time independent RTF as being convenient. Consequently, the 
time dependence of the radiation field originates exclusively from the coupling to the FLASH code using the energy 
source term from Fquation 

In this section, we show results from testing the evolution of the source term by following the propagation of the radi¬ 
ation field in a highly opaque medium. The source function is updated using a simple forward Fuler time integration 
of the energy source term. Since the radiation field shows a diffusion like evolution in the limit of high optical depths, 
we compare our numerical solution to the analytic solution of the diffusion equation. 


3.3.1. Setup 

In this test, we investigate the ability of our solver to follow the flux of radiative energy into a highly opaque 
medium. In this case, the propagation of the radiation field can be described by the diffusion approximation, and we 
show that our approach reproduces the diffusion limit accurately. The diffusion approximation is derived from the mo¬ 
ment equations of the RTF by invoking a closure relation between the radiative energy and the radiative pressure (e.g., 
the Fddington approximation). The radiation equations themselves then form a hyperbolic system. By neglecting 
the explicit time dependence of the radiative flux F and assuming that F oc VF^, the flux can be eliminated from the 
equations . The dynamics of the radiation fiel d J - cEffiAn) can then be described in a single equation, the diffusion 
equation ( Mihalas and Weibel Mihalasl 1984 ): 




dt 


\nx 


SJJ\^cx{S-J). 


(30) 
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Figure 7: The solutions of the IPascucci et al. benchmark problem. Left column: the temperature structure through the 

xz-midplane of the disc for total radial optical depths of t = 1 (top), t = 10 (mid), and t = 100 (bottom). Right column: averaged 
temperature profiles in the xy-midplane in comparison with the solutions of RADMC-3D. Solutions obtained with FLASH/RT 
use 768 directions for the angular discretization. Monte Carlo computations with RADMC-3D were performed using 10* photon 
packages. 
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where n denotes the number of dimension. We do not allow any interaction of the radiation field with the hydrody¬ 
namics and only follow the propagation of the radiation field. Hence, the diffusion equation becomes homogeneous 
since S = 7. In this case, the solution to the diffusion equation is described by the Gaussian function 


7/5 (x, t) = 


Jq 




exp 


(X - Xq)^ 

ADt 


(31) 


where 7o denotes the initial mean intensity at t = to, xq is its initial position, and D - cKrix) is the diffusion 
coefficient. We use Equation (ISTT i to compute the initial conditions 7(x, to) for our test setup. We perform ID and 3D 
computations with the initial conditions Jo = 7(xo, to) - 10^ erg s“' cm“^ sr“* with to - 10“'* s in 3D and to = lO^*" s 
in ID respectively. The center of the Gaussian is at xq = 0, and we evolve the radiation field until t = 20 x fo is 
reached. The length of the computational domain is 1 cm with a homogeneous density distribution of p = 1 g cm * and 
a constant absorption coefficient k = 1000 cm^ g“*, which results in a highly optically thick medium. The temperature 
is constant and arbitrarily set to T = IK. Since no heating or cooling is allowed, there is no hydrodynamical response 
from the medium and all hydrodynamical quantities are constant in space and time. 

Since we solve the time-independent RTE, there is a problem in reproducing the time-dependent term in Equation 
(l30l l. Strictly speaking, the static source function vanishes since we do not couple the radiation field to the medium 
through which it propagates. Consequently, the mean intensity would also vanish in the ti me independent so lution. 


However, the time dependence causes an effective contribution in the source function (e.g. iJack et al.L 120121) if the 


time discretization is carried out implicitly in the RTE ([TJ. This contribution depends on the specific intensities of the 
previous time step, is evolved through time, and describes the evolution of the radiation field. Since we do not account 
for this implicit contribution (which would require to store the complete scalar field of angle dependent specific 
intensities), we solve the problem by operator splitting using the right-hand side of Equation (l30l l to calculate the 
new source function at the following time step. The evolution is done using a simple forward Euler time integration 
scheme of the form 

Sn ^ S„-i + At„xc(J„-i - Sn-l), (32) 

where Af„ is the length of the current time step n. Therefore, the time step is restricted to be (Section|Z4Tj 


Af„ = min(--j ^ j ( 33 ) 

where the min function denotes the minimum change in the source function with time from all cells i in the computa¬ 
tional domain (ELASH does not support adaptive time stepping on a block level but rather uses a uniform global time 
step), krad limits the maximum change in the source function, and we found a value of krad ~ 0.1 to give stable and 
accurate results in 3D. 


3.3.2. Results 

The results of the ID solutions are shown in Eigure[8] We compare the numerical results with the analytic solution 
given by Equation (OTT i and found our results to be within 1 % accuracy at a resolution larger than 32 cells. At the edge 
of the domain, the numerical solution deviates from the diffusion solution as radiative energy can leave the domain and 
we allow no irradiation from the outside. The results from the 3D computation are shown in Eigurej^and compared 
to the diffusion solution along the three main axes of the domain. In the 3D case, the domain is subdivided by the 
AMR grid into 4 blocks in each dimension. Each block contains 8^ cells giving a total linear resolution of 32 cells. In 
the 3D case, the setup consists of a Gaussian kernel around the origin which diffuses outwards. The solutions along 
each coordinate axis are obviously indistinguishable, emphasizing the accuracy and importance of the homogeneous 
angular HEALPix tessellation (AppendixB|i. The 3D computations were performed using 192 directions. 


3.4. ID Non-equilibrium Radiative Shocks 

Testing the radiative transfer solver for radiative shock computations is the next crucial step and requires the 
combination of our radiation solver with the ELASH code. The source term is determined by the energy budget of 
absorption and emission processes. We recall the frequency integrated source term from Equation ([8]l here: 


Grad = ‘^ttXaiJ “ B), 
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Figure 8: Results of the ID diffusion test for different homogeneous spatial resolutions, top: nx=16, mid: nx=32, right: nx=64. 
The dashed lines show the initial conditions at f = to determined by the Gaussian solution of the diffusion equafion. The initial 
radiative energy (symbols) is evolved and diffuses outwards until f = 20 x to is reached and compared to the analytical solution 
(solid lines) of the homogeneous diffusion equation. For a sufficient spatial resolution, the numerical solution stays within 1% 
accuracy. At the edge of the domain, the radiation solver deviates from the diffusion solution as radiation leaves the domain while 
the diffusion solution is valid for an infinite domain. 






error along y-axis 






x[cm] 


y[cm] 


z[cm] 


Figure 9: Results of the 3D diffusion test along the x-(left), y-(mid) and z-axis (right) of the simulation box with a homogeneous 
spatial resolution of nx=ny=nz=32 (symbols). Top row: The dashed lines show the initial conditions at f = fo determined by the 
Gaussian solution of the diffusion equation. Blue solid lines show the analytic solution according to Equation[3T] Symbols show 
the numerical results from our radiation solver. Bottom row: The relative error of the numerical solution. The 3D solution is not as 
accurate as the ID results but still within 10% of the analytical solution. The obvious independence of the solution on the direction 
axis results from the homogeneous angular HEALPix tessellation. The calculations were performed using 192 directions. 
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which is coupled to the hydrodynamical solver by adding it to the right-hand side of the Euler equation for the internal 
gas energy. For this test case, the emission and absorption opacities are equal. Since the shock setup is used for test 
purposes, we neglect the magnetic field. 


3.4.1. Initial Conditions 

The initial conditions are consistent with the theoretical work of iLowrie and EdwardsI (l2008h . In their work, the 
jump conditions and the equations of radiation hydrodynamics are given in a dimensionless form. The equations are 
normalized using reference material quantities and a constant Pq which arises from the normalization process and is 
given by 


Po = 


ml' 


(35) 


The quantities denoted with a tilde are the dimensional reference material attributes (temperature Tq, density po, sound 
speed So) and Oy is the radiation constant. The ”0”-subscript indicates pre-shock state initial values. Pq gives a measure 
for the relative importance of gas an d radiation pressure or alternatively, the radiative energy to the material energy 
(IMihalas and Weibel Mihalasill984l) . For our test setups, we choose a grey non-equilibrium shock setup with Mach 
numbers of Mq = 1.2, Mq - 2 (subcrit ical), and Mn = 5 (supercrit ical), which we compute in the reference frame of 
the shock with Pq = 10 and 7 = 5/3. Lowrie and Edwarda l 2008h give a dimensionless absorption and transmission 
cross section, which determine the radiative energy exchange and diffusivity of the radiating materials. Evaluating the 
dimensionless values gives an absorption coefficient of A'fl » 423.0cm^/g and atotalextinctioncoefficientof;i^ 788.0 
cm^/g, which results in an effective photon destruction probability of e = Kalx ~ 0.5377. The initial dimensionless 
pre-shock gas temperature Tq and density po are set to unity, the post-shock initial values (Ti, pi) are computed using 
the Rankine-Hugoniot jump conditions. The actual dimensional initial conditions can t hen b e calculated using their 
dimensional reference material values (for more details we refer to lLowrie and EdwardsI (l2008h '). Finally, the radiation 
temperature 


7 ’,.= 


J 


C^SB 


1/4 


(36) 


is initially in equilibrium with the gas temperature. For the radiation shock test problem, the source function is 
determined by a thermal emission and a diffusive part. This is equivalent to using the isotropic scattering source 
function 

S^{l-e)J + eB (37) 


with the appropriate photon destruction probability and a thermal energy contribution given by the frequency inte¬ 
grated Planck emission B - p4 sijjce the radiation field will not be not be in thermal equilibrium with the 
material throughout the simulation, we need to iterate until a consistent solution of the mean intensity J is found. 
However, since e 0.5377 gives only a moderate scattering contribution and using the solution from the previous 
time step, the accelerated lambda iteration usually converges after 2 or 3 iteration steps. 


3.4.2. Results 

The shocks need a few nanoseconds to relax into a static equilibrium state. Figure[T0l[TT]and[T2]show the resulting 
temperature and density profiles after 10 nanoseconds. Sufficiently far upstream (left) and downstream (right) of the 
hydrodynamical shock (at x - 0), gas and radiation are in a thermodynamical equilibrium and the radiation tempera¬ 
ture coincides with the gas temperature computed from the initial conditions. Since the total extinction coefficient 
is about twice the thermal absorption and emission coefficient, the temperature of the radiation field and the gas are 
out of equilibrium near the shock front. 

The subcritical shock with Mq = 1.2 (Figure [TOli shows a hydro shock but no spike in the radiation temperature. For 
Mq = 2 the so called Zel’Dovich spike in the gas temperature appears for the first time as seen in Figure fTTI The spike 
appears since radiation is transported through the hydrodynamical shock and preheats the inflowing gas, which is 
initially in a thermal equilibrium with the radiation field in the upstream region. After the gas has passed the hydrody¬ 
namical shock, it cools down until the radiation field and the gas are again in thermal equilibrium on the downstream 
side of the shock. Since the upstream temperature at the shock front is still less than the downstream temperature 
the shock is subcritical. For Mq = 5, the shock becomes supercritical, since the upstream gas is preheated until it 
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Figure 10: Normalized temperature and density profiles for the subcritical shock with Mo = 1.2 in the equilibrium state after 10 
nanoseconds. The gas is preheated on the upstream side and cools downe on the downstream side of the hydro shock front. 
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Figure 11: Same conditions as in Figure [TO] but with Mo = 2. The maximum temperature at the shock begins to exceed the 
downstream equilibrium temperature which results in the Zel’dovich spike. Since the temperature at the upstream side of the shock 
is still well below the downstream temperature, the shock is subcritical. 


reaches the downstream gas temperature even before passing the hydrodynamical shock front. The discontinuity in 
the gas temperature is then restr icted to the narrow range of the Zel’dovich spike (Figure [TSTi. Our solutions resemble 
the semi-analytical results from lLowrie and Edwards! (120081) and show the correct spike evolution. However, a closer 
look at the results show a slight deviation of the shock front from its initial position (at x - 0). Especially in the 
supercritical case, the shock front drifts very slowly into the downstream direction. This drift is due to the absence 
of the radiation pressure in our approach, which becomes important for high Mach numbers (with a high downstream 
gas temperature). While the shock front drifts very slowly, the temperature and density profiles do not change since 
the radiation source term is still very well approximated in our approach. 


4. 3D Collapse Simulations 

In this section, we show results from full 3D radiation hydrodynamical simulations performed with FLASH/RT. 
Since we aim to use our framework for the modelling of radiative feedback in star formation simulations, we show the 
capabilities of our method in two self-gravitating collapsing cloud simulations. We follow the collapse until the first 
hydrostatic core is formed and before the dissociation of hydrogen molecules start (the first collapse). In Section |472] 
we show results from a basic collapse simulation without rotation and compare the resulting profiles to other similar 
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Figure 12: Same conditions as in Figure[TO]but with Mq = 5. The temperature on the upstream side of the hydro shock front 
reaches the downstream equilibrium value. The Zel’dovich spike gets very narrow and the shock becomes supercritical. 


works. Afterwards, we show results from a more complex simulation including rotation and turbulence (Section l43T l 
and compare the results to a simulation without modelling radiative transfer. The angular resolution of the radiative 
transfer calculations are the same for both collapse simulations, and we use 768 directions to compute the radiation 
held (nSide=8 for the HEALPix tessellation). 


4.1. Opacities 

Since our solver does not yet support any frequency dependence, the source function S is only determined by 

CT T'^ 

the frequency-integrated thermal emission of the gas (S — B — —), and we neglect any scattering processes. 
Consequentl y, we have to us e frequ ency-integrated mean dust opacities. For this purpose, we choose the Planck mean 
opacities bv ISemenov et al.l (l2003h . In their work, the dust composition model takes into account the evaporation 
temperatures of ice, silicates, iron as well as their density dependencies. We coupled their subroutine^ for computing 
temperature and density dependent dust opacities into FLASH, and we choose the input parameters for spherical 
homogeneous dust grains with a normal relative iron content in the silicates of Fe/(FeH-Mg) = 0.3. 


4.2. Collapse without Rotation 

In this section, we study the collapse of a spherical, homogeneous, and gravitationally unstable density distribu¬ 
tion. The initial conditions do not contain any turbulence or density perturbations and hence, the results are spherically 
symmetric. This setup represents a common benchmark for the capabilitie s of a radiation hydrodyn a mical astrophys- 
ical co mputer code, and we compare our results t o sim ilar work done by Commercon et al. (2011), Masunaga et al.l 
(Il998h . and the pioneering simulations of lLarsonl(ll969h . 


4.2.1. Initial Conditions 

We start with highly gravitationally unstable initial conditions. The cloud core of one solar mass consists of a 
homogeneous sphere with radius /?o = 7.07 x 10'® cm (a; 4725 AU) and and density po = 1.38 x 10“** gcm“^, which 
results in an initial free fall time of fff a; 56.67 kyrs. The linear size of the 3D computational domain is four times the 
initial cloud radius Rq in each dimension. The surrounding gas density is a hundred times less than the initial cloud 
density po, and the cloud is initially in thermal equilibrium with the ambient gas at a temperature of To = 10 K resulting 
in an initial isothermal sound speed of Cs ~ 0.195 kms“’. Since the cloud is initially not in pressure equilibrium with 
its surroundings, FLASH’S hydrodynamical solver drives a weak shock wave into the ambient gas which is soon 
dissipated. To prevent our radiation solver from resolving this shock in terms of radiative energy exchange (which 


^http://www.mpia-hd.mpg.de/homes/henning/ 
Dust.opacities/Opacities/opacities.html 
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Reference 

Rfc [AU] 

Mfc [Mol 

Tfc[K] 

Tc[K] 

This work 

4 

1 X 10-^ 

50 

186 

Commercon et al. 120111 

8 

2.1 X 10-2 

81 

396 

Masunaga et al. 119981 

8 

^ 10-2 

60 

200 

Larson 119691 

4 

1 X 10-2 

- 

170 


Table 1; Comparison of simulation results; Rfc is the radius of the first core, is the core mass, Tfc the central 
temperature and is the temperature at R^. 


would result in rather small time steps), we do not couple the radiation field to the hydrodynamics outside of Rq but 
rather keep the ambient gas and radiation temperature fixed. 

The initial conditions result in a gravitationally unstable cloud core which contains nearly two Jeans ma sses. To ensure 


a prop er resolution and avoid artificial fragmentation during the collapse, we use the Jeans condition bv lTruelove et al 


(Il997 ) as the refinement criterion of the AMR grid. In our case, we use at least Nj = 9 grid cells per Jeans length. To 
resolve the first hydrostatic core properly, we allow a maximum linear resolution of Ax « 0.07 AU which requires the 
AMR grid to cover 11 levels of resolution. 

The summarized initial conditions are; 


Mass 

M 

Density 

Pit 

Temperature 

To 

Angular Velocity 

Q 

Radius 

Ro 

Free Fall Time 



1.0 Mo, 

1.38X 10^'*gcm^^ 
10 K, 

O.Orads ', 

7.07 X 10'® cm, 
56.67 kyrs. 


4.2.2. Results 

The cloud core starts to collapse, and as soon as the maximum density in the cloud exceeds about 10 gcm“^, 
the central regions of the cloud core become optically thick. At this point, the central temperature starts to rise rapidly 
and the following evolution proceeds almost adiabatically with more gas falling onto the central quasi-hydrostatic 
core. Since the simulation does not contain any rotation or turbulence, the 3D solution is spherically symmetric, and 
we present the results in the form of averaged radial profiles. The profiles for density, radial velocity, temperature, 
optical depth, and central mass after 1.036 x ts are shown in Figure [T3] The resulting protostellar core has a mass of 
Mfc 1 X 10“^ Mq, a radius of ~ 4 AU, and a central temperature of Tc » 186 K. The boundary of the core can be 
identified easily in the velocity profile, where there is a sudden decrease in the infall velocity (the accretion shock). 
Inside the core, the infall does not stop completely in dicating that th e core is only quasi-hydrostatic. 

Our result s are quantitatively very simi lar to those of ( j_969l) and qualitatively very similar to the more recent 

works by Masunaga et al. (1998) and Commercon et al. (I 2 OI 1 ). Table [1] shows an overview of the characteristic 
temperature, mass and radius of the first core in comparison to these works (the common reference point is when the 
maximum central density of the first core reaches pfc ~ 2 x lO^’^gcm^^). Apparently, our computations produce 
qualitatively similar results, although the methods invoked in the other works are quite different and use different 
initial conditions and opacity models. 


4.3. Collapse with Rotation and Turbulence 

This simulation run has very similar initial conditions as described in the previous section except that we add 
rotational and turbulent energy. The cloud is initially in a rigid body rotation around the z-axis at the center of the 
simulation box. The ratio of rotational and gravitational energy is given by 


1 /?3 q2 

3 GMo' 


(38) 
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Figure 13: Profiles of the collapse simulation after t = 1.036 tg, the maximum density at the core center is pc ~ 2.0 x 10 
with a temperature of Tc ~ 186 K, a radius of R(c ~ 4 AU and a mass of Me ~ 10“^ Mo 








We choose jS = 0.03 which gives an initial angular velocity of Q n = 1.886 x 10 '^rads ' and agrees with typically 


observed values of molecular cloud cores (iGoodman et al.L [19931). In addition, we superimpose a turbulent velocity 


perturbation on the initial uniform angu lar velocity field. Th e construction of the velocity perturbation is based on the 
theory for incompressible turbulence bv lKolmogorovl(ll94ll) . in which the kinetic energy E of the velocity fluctuation 
with wave number k is described by a power spectrum 


E{k) cx kP. 


(39) 


The wave number k — Injl is the inverse of the length scale / of a turbulent fluctuation (sometimes called eddy). In 
our case, the spectrum has a power law index of p = -2 resembling a Burgers type model of turbulent energy decay. 
The geometries and density distribution of the initial cloud core are the same as for the simulation without rotation 
and turbulence. 

In addition to the simulation run with FLASH/RT, we also run the simulation without modelling radiative transfer. 
Instead, we use a barotropic EOS with a density-dependent effective adiabatic exponent y that mimics radiative cool¬ 
ing. The internal energy/temperature is fixed at To = 10 K as long as the gas density is less than p « 10“^^ gcm“^ 
(isothermal). Above this threshold density, the temperature rises slowly with y - 1.1 until the adiabatic exponent 
becomes y - 4/3 above p ^ 10^'^gcm^^ (adiabatic). We ran the simulation including radiative transfer as well 
as the reference run with the barotropic EOS until the formation of the first hydrostatic core with a central density 
of pfc ~ 10 " gcm~^. At this point, both simulations cover 9 different levels of resolution in the AMR grid with a 
maximum linear resolution of Ax » 0.57 AU while the whole simulation box has a linear extent of 18903 AU. 

The summarized initial conditions are; 


Mass 

M = 1.0 Mo, 

Density 

po = 1.38 X 10“^* gcm“^, 

Temperature 

To = 10 K, 

Angular Velocity 

Q = 1.886 X 10^^^ rad s^ 

Rotational Energy 

P = 0.03, 

R = 7.07x 10'® cm. 

Gravitational Energy 

Radius 

Eree Eall Time 

fff = 56.67 kyrs. 


4.3.1. Results 

The rotational energy and the superimposed turbulent velocity perturbations break the symmetry of the simulation. 
Eigure[T4| shows the column densities along the main axes of the inner region where the dense first core has formed 
after about 60kyrs (a; 1.07 fff) with a maximum gas density of a; 10“'^ gcm“^. Because of the additional rota¬ 
tional and turbulent energy, the formation of the first core is deferred and forms later than in the previous simulation 
(Section [4.2l) . The conservation of angular momentum causes the first core to be flattened roughly along the z-axis 
and the density distribution shows a flat disc-like structure revolving around the central compact hydrostatic core. 
The resulting density distribution is roughly the same as in the reference run without radiative transfer. The initial 
collapse which seeds the formation of the central core does mostly occur in the isothermal phase, hence, modelling 
radiative feedback does not influence the initial formation of the core significantly. However, Eigure [15] shows the 
resulting density weighted temperature averages along the main axes in the central regions around the first core (e.g 
Jp r dz! f pdz). The left column shows the results including radiative transfer (ELASH/RT) while the right column 
shows results from the reference run. The ELASH/RT model clearly shows h ow the central core hea ts the surround¬ 
ing gas to a temperature roughly 30% higher than in the reference run (like in IPrice and Bate! ( 2010l) '). The resulting 
temperature density distribution in comparison to the barotropic EOS is shown in Eigure [Thl 

Unfortunately, our ELASH/RT simulations are very costly (see Section[5]for more details) and currently, it is not fea¬ 
sible to continue these simulations with out coupling the radiative transfer solver to a subgrid model for the formation 


of the central core, e.g., sink particles ( Eederrath et al. . 2010l) . However, our current test simulations show the first 


stages of disc formation and the importance of modelling radiative transfer accurately. Since the thermodynamics of 
the gas significantly influence the fragmentation behaviour, modelling radiative transfer is indispensable to study the 
further evolution of the protostar, the circumstellar disc, and the surrounding gas envelope. 
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Figure 14: Column densities along the z- (left) and y-axis (right) of the simulation box after the formation of the first hydrostatic 
core at r « 60 kyrs « 1.07 tff including rotation and turbulence. The rotational energy forces the gas to accumulate in a circumstellar 
disc (in the xy-plane) around the first core. 
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Figure 15: The plots show density weighted temperature averages (e.g., f pT dz! f pdz) from a collapse calculation including 
rotation and turbulence. Left: Results from the FLASH/RT calculations including radiative transfer. Right: Results from FLASFl 
calculations using a barotropic EOS. The ambient gas temperature in the FLASH/RT models is about 30% higher. 
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nr. of cores 

Time [s] 

Speedup 

Blocks per cpu 

Performace per Block [%] 

96 

86.06 

1.0 

37-39 

100.0 

144 

60.60 

1.42 

25-26 

95.2 

192 

48.01 

1.79 

18-20 

89.6 

240 

41.11 

2.09 

14-16 

82.6 

288 

34.93 

2.46 

12-13 

81.0 

336 

32.98 

2.60 

10-12 

75.5 


Table 2: Results from the scaling test normalized to a run with 96 cpus; because of the increased communication overhead, each 
cpu should handle as many block as possible in terms of memory requirements. 


5. Performance 


The FLASH code shows excellent scaling behaviour on any computational infrastructure (e.g. Frvxell et al. , 2000l) . 


For this work, the computations are clearly dominated by the solution of the RTF. Hence, the scaling behaviour of 
the radiative transfer solver is crucial for the total performance of the FLASH/RT calculations. We investigate the 
scaling performance of our radiation code using the disc benchmark setup (see Section lT^ . We performed 50 formal 
solutions of the RTF using 192 directions on a spatial range covering 5 refinement levels. After the initial refinement, 
depending on the density structure and the radius, the computational domain consists of 3648 valid subdomains (leaf 
blocks) each containing 8^ cells. The FLASH code distributes the blocks among all available MPI ranks using a 
Morton space-filling curvfl The scaling tests were run at The North-German Supercomputing Alliance in Berlin on 
the Cray XC30 ’’Gottfried” using 12-core Xeon IvyBridge processors. Figure[T2]and Table|2]show the scaling results 
for the computation of the formal solution of the RTF averaged over 50 cycles. The scaling is normalized to the wall- 
clock time using 96 cores (e.g, 8 Xeon IvyBridge processors). ’’Gottfried” provides 2 Xeon processors with 24 cores 
in total per computing node, hence, adding 24 cores to the computation will increase the communication overhead. 
Figure fTTbl shows the speedup compared to a perfect scaling behaviour. The radiation solver scales reasonably well 
considering the communication of non-local information, which is necessary for the solution of the RTF. Figure [TTcl 
shows that doubling the number of cores decreases the performance per block by approximately 10%, which we 
consider also as reasonable. 

The cost of the radiative transfer solver from a 3D collapse simulation (Figure flTdl l is comparable to the cost for the 
computation of the self-gravitational potential which is done by a Poisson tree-solver. However, the radiative transfer 
solver in this par ticular simulation u ses a rather moderate angular resolution of 192 directions (using the HFALPix 
tessellation from Gorski et al. ( 2005h '). For runs including rotation and turbulence, the angular resolution probably 
needs a much higher resolution of at least 768 directions or higher. Since the cost of the radiative transfer solver 
scales linearly with the number of directions, it dominates the entire simulation run compared to the calculation of 
self-gravity. So far, we have tested the FLASH/RT code on our own computing cluster in Hamburg (32 nodes with 
2x Intel Xeon Hexa-Core CPUs, 2.40 GHz) and at the North-German Supercomputing Alliance in Berlin on the Cray 
XC30. 


6. Summary 

We have implemented a new radiation transfer solver based on the method of hybrid characteristics. The solver 
successfully reproduces standard radiative transfer problems, including NLTF, thermal radiative transfer and the dif¬ 
fusion limit. We proved the feasibility of the method for 3D collapse simulations where radiative transfer is the 
dominant cooling process during the formation of the first protostellar core. In contrast to the FLD approximation, 
our method preserves the anisotropy of the radiation field which becomes crucial in the transition from optically thin 
to optically thick regions (e.g, a protostellar disc). The radiation solver is implemented in the framework of the MHD 
code FLASH which allows for a straight forward coupling of both codes (e.g., the collapse simulations. 

However, the explicit energy coupling, as described in Section l2.4.11 puts rather strong limitations on the time step. 


^http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug/ 
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A possible improvement can be achieved by combining the raytracer with the solution of the moment equation for 
the radiative energy. In contrast to the FLD appr oach, one can com pute an angle dependent diffusion coefficient in 


the form of the variable Eddington tensor (VET) (I.Tiang et al.L 1201 2h which can be achieved using our raytracer. The 


advantage is that the evolution of the radiati ve energy can be handled implicitly by solving the linearized moment 
equation for the radiation temperature, like in ICommercon et al.l (1201 ih . which res olves the probl e m of t he time step 
restriction. The framework for this has already been implemented in ELASH by iKlassen et al.l (120141) and can be 
combined with our raytracer to implement the VET approach. 


Our m etho d is a generalized a nd enhanced implementation of the hybrid characteristics method by iRiikhorst et al 


( 20061) and Peters et al. ( 2010l) . The original implementation was restricted to direct irradiation from point sources 
and the integration of optical depths respectively. The ELASH code in combination with our radiative transfer frame¬ 
work allows for the solution of a much wider range of problems and can also very easily be extended to handle a 
more complex form of the radiative transfer equation. Our implementation fits very well into the parallel design of 
the ELASH code which is based on AMR with domain decomposition. Our method works within the AMR design of 
ELASH and is able to solve the 3D RTE on a wide range of scales which is indispensable for star formation simula¬ 
tions. 


Acknowledgments 

L.B. acknowledges financial support by the Deutsche Eorschungsgemeinschaft (DEG) mainly via the Emmy 
Noether grant BA 3706/1-1 Theory of Massive Star Formation and partly through the Graduiertenkolleg 1351 Extraso¬ 
lar Planets and their Host Stars and the Priority Program 1573 Physics of the Interstellar Medium. T.P. acknowledges 
financial support through a Eorschungskredit of the University of Zurich, grant no. PK-13-112, and from the DEG 
Priority Program 1573 Physics of the Interstellar Medium. This work benefited from helpful discussions with Peter 
Hauschildt (Universitat Hamburg) and Stefan Dreizler (Universitat Gottingen). Most of the collapse simulations were 
carried out at The North-German Supercomputing Alliance in Berlin {Gottfried, HERN). 

AppendixA. Accelerated Lambda Iteration 

The lambda operator A describes the task to compute the radiation field from the source function. It is usually 
written as 

7 = A[5]. (A.l) 

Eormally, we can solve this by inverting the Lambda operator. When we arrange the cells of a 3D domain successively 
in a ID vector, we can write the operator as a matrix. But the complete operator for one cell in the computational do¬ 
main contains all radiative contributions from each other cell. Hence, the Lambda matrix is far from being sparse. The 
explicit construction and storage of the Lambda matrix would easily reach computational limits in terms of memory 
requirements. Eurthermore, the inversion of the Lambda operator is far too costly to be used in 3D radiative transfer. 
Instead, the formal solution (|4| is used. Since the source function may depend on the mean intensity, this task requires 
iteration over Equations (|2l) to (|4](. This is called Lambda iteration but it usually fails in optically thick regimes. This 
happens because photons can be trapped and scattered many times, if a single cell of the computational domain is 
optically thick. The ordinary Lambda iteration is not able to account for these processes on scales smaller than the 
spatial resolution. 

The idea behind the accelerated Lambda iteration (ALI), is to extract these sub-cell scattering contributions from 
the Lambda operator (and hence from the iteration), because we are not able to resolve them anyway. The ex¬ 
tracted part of the ordinary Lambda operator is then put into a new approximated Lambda operator, which is solved 
quasi-analytically. Since the approximated Lambda operator usually only contains a small part of the whole Lambda 
operator (the subgrid part so to say), it is easy to compute, store and solve. Mathematically, the Lambda operator 
becomes split 

A = (A - A*) -H A* (A.2) 
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where A* denotes the approximated Lambda operator. Inserting this into eauation lA.il and using the source function 
for isotropic scattering (Eauationl22li. we get 

S ^eB + (l- e)(A - A*)5 + (1 - e)A*5 (A.3) 

Since the A*-operator consists of only a small part of the whole Lambda-operator, it is sparse and easy to solve. We 
bring it to the left-hand side: 

[1 - (1 - e)A*]5 ^eB + (l- e)(A - A*)S. (A.4) 

We introduce the iteration scheme, because there is still a contribution of the source function on the right-hand side. 
This remaining contribution can be regarded as the non-local contribution of the radiation field, which is solved by 
iteration. Inverting the approximated Lambda-operator then yields 

= [1 - (1 - e)A*]-\eB H- (1 - e)(A - A*)5"). (A.5) 

The scheme in Equation |A3] is a combination of iteration and analytic solution. The non-local contributions (in the 
Lambda matrix (A - A*)) are accounted for by iteration while the local subgrid scattering is handled by an inversion 
of the approximated Lambda operator (A*). The computational cost of the inversion of the A*-operator depends on its 
bandwidth, which determines the range on which we solve analytically. Obviously, a diagonal A* is trivial to invert. 
But since the diagonal part of the Lambda operator describes only the local scattering in a single cell, it is not the 
best choice in terms of iterative performance. Usually, a tri-diagonal operator yields the best compromise between 
fast convergence and computational cost. But this requires the solution of a coupled set of linear equations, which is 
complex to implement. Eor now, we stay with a diagonal local A*-operator, since it is the easiest one to implement 
and still has a tremendous effect on the convergence rate. 


AppendixB. The Angular Discretization using HEALPix 


The choice of the solid angle grid is equivalent to the problem of discretizing the surface of a unit sphere. The 
method of characteristics requires the solution of the parameterized RTE along a large number of directions n de¬ 
pending on the anisotropy of the specific intensity /(x, n). In general, this requires a homogeneous dis cretization of 
the so lid angle Q on the 4n unit sphere. Eor this purpose, we use the HEALPijJI scheme introduced by Gorski et al 


(l2005l) . HEALPix ensures an optimal discretization of the unit sphere (also called pixelation or tessellation) into a 
number of finite solid angles AQ. HEALPix in general addresses problems in which a function on domains of spher¬ 
ical topology has to be analyzed. The pixelation scheme was originally developed to handle large datasets generated 
by cosmic microwave background experiments (e.g., WMAP, Planck) and provides a software librar}@ with numerous 
subroutines for spherical discretization and numerical analysis of functions or datasets on the sphere. 

The HEALPix pixelation has a base resolution of 12 pixels in three rings around the poles and the equator of the unit 
sphere each covering the same area. Based on these base pixels, the resolution is refined by dividing each base pixel 
into subpixels, where Aside has to be a power of 2 (Aside = 1,2,4,8,...). The total number of pixels (assuming an 
isotropic refinement) is then Apix = 12A^jjg). 


^Hierarchical Equal Area isoLatitude Pixelization 
^ http://healpix.j pi.nasa.gov/ 
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Figure 16: Temperature distribution with respect to the gas density in the simulation box at the end of the collapse simulation 
including rotation and turbulence. Black dots show the temperature distribution from the FLASff/RT run, red dots resemble the 
temperature density dependence of the barotropic EOS. 
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(a) Wallclock Time 
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Figure 17: Results from the parallel scaling test; In a) the total wall-clock time for the formal solution averaged over 50 iteration 
cycles from the Pascucci disc benchmark is shown. In h) the speedup normalized to the wall-clock time using 96 cores is shown. 
The number in brackets denote the minimum and maximum blocks per core, which is the result from the Morton space-filling 
curve. In c) the performance per block is shown, which decreases by roughly 10% if the number of cores is doubled. In d) the 
fractional runtimes for the most costly steps of the the collapse simulation (Section|4j2} is shown. All computations were performed 
using 192 directions. 
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Figure B.18: The FIEALPix tessellation scheme, from iGorski et alj ( l2005h . The Area in light grey shows one of the eight (four 
north, and four south) polar base pixels and the dark grey area shows one of the four equatorial base pixels. Moving clockwise 
from the upper left panel the base pixels are hierarchically subdivided with the grid resolution parameter equal to Aside = 1, 2 , 4, 8 
and the total number of pixels Apix = 12,48, 192,768. 
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